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We illustrate the stochastic method for solving the Schwinger-Dyson equations in large- A'^ quantum 
field theories described in ArXiv: 1009.4033 on the example of the Gross- Witten unitary matrix 
model. In the strong-coupling hmit, this method can be applied directly, while in the weak-coupling 
limit we change the variables from compact to noncompact ones in order to cast the Schwinger-Dyson 
equations in the stochastic form. This leads to a new action with an infinite number of higher-order 
interaction terms. Nevertheless, such an action can be efficiently handled. This suggests the way 
to apply the method of ArXiv: 1009.4033 to field theories with U (N) field variables as well as to 
effective field theories in the large- A'^ limit. 

PACS numbers: 02.70.-c; 02.50.Ey; 11.15.Pg 

Keywords: simulation methods, planar approximation, lattice field theory, random processes 



Introduction 
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In many physically interesting situations direct Monte- 
Carlo samphng of the path integral of the theory is not an 
efficient simulation method. Probably the most notable 
example in hadronic physics is QCD at finite chemical 
potential. Another example are quantum field theories 
in the large- limit, since the computational complexity 
of Monte-Carlo algorithms tends to infinity in this limit, 
and one always has to rely on extrapolations from finite 
N. 

One of the ways to overcome these difficulties is to ex- 
pand the desired path integral in powers of some coupling 
constants, and to sum up the obtained series stochasti- 
cally. This is the basic idea of the Diagrammatic Monte- 
Carlo [H, This method, however, cannot be straight- 
forwardly applied to field theories with continuous vari- 
ables, since in such cases the perturbative series are typ- 
ically only asymptotic Q. 

For field theories with a non-Abelian internal symme- 
try group, such as U (N), each diagram contains also 
some power of N, and the sums over the terms leading 
in TV are typically convergent Q. This salient feature of 
large- iV quantum field theories was used in [5J to con- 
struct a stochastic algorithm for summing over all planar 
diagrams in perturbative expansions. This algorithm is 
based on the stochastic solution of the Schwinger-Dyson 
(SD) equations in the factorized form, so that they are in- 
terpreted as the equations for the stationary probability 
distribution of a certain random process. The correlators 
of the field variables are thus estimated as probability 
distributions of some random variables, similarly to the 
"worm algorithm" [6]. Once the SD equations can be 
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where w {x) schematically denote the field correlators and 
the coefficients Pc (x), Pe ix\y) andpj 2/2) satisfy the 

inequality be {x) \ + \pe {x\yi) \ + \pj {x\yi,y2) \ < 1, one 

can construct a "nonlinear" random process for which 
w {x) is the stationary probability distribution Q . 

However, when the field variables are the elements of 
U {N), SD equations can be rewritten as ([T]) only in the 
strong coupling limit. On the other hand, the continuum 
limit of such theories typically corresponds to the weak- 
coupling limit. In this paper we consider the simplest 
nontrivial model defined by the integral over U (N) group 
in the limit N ^ 00, namely, the Gross- Witten unitary 
matrix model Q- We demonstrate that the method of 
[1] can be used to simulate this model in both the strong- 
and the weak-coupling regimes. In the latter case, one 
has to change the variables to the noncompact ones, 
which are more adequate for the weak-coupling expan- 
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FIG. 1; Observables Gn (A) in the strong-coupling regime. 



Strong-coupling regime 

The SD equations for the model ^ in terms of the 
observables G„ (A) are: 

Gi + iG2-i-0(3) 
G„ + ^ Gfc Gn-k + J G„+i - - G„_i =0, n > 2 (4) 
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By a variable redefinition G„ = J\fc"wn one can trans- 
form these equations into the stochastic form ([1]) if A > 
As « 5.94. According to [Hj, such equations describe the 
stationary distribution of the random process operating 
on the stack of the elements of the form {n, ±}, where 
n > 1 is an integer labeling the observables G„. Wn is 
proportional to the difference of the probabilities to find 
the elements {n, -I-} and {n, — } at the top of the stack. 
At each step of the random process one performs at ran- 
dom one of the following steps: 

Create : With probability {X\fc)~^ push a new element 
{1, -I-} to the stack. 

Join : With probability TV pop the two elements 
{m, si} and {71,82} and push a new element 
{m + n, -8182}- 

Increase : With probability (Ac)~^ increase the top- 
most element in the stack by one. 

Decrease : With probability c/A decrease the topmost 
element in the stack by one and change its sign. 

Restart : Otherwise restart with a stack containing a 
single element {1,-|-}. 

The dependence of the observables G„ (A) on A is il- 
lustrated on Fig. [T] The data were obtained for 10^ 
iterations of the random process, which took several sec- 
onds on a standard 2 GHz CPU. In agreement with the 
exact solution 7], only Gi (A) is nonzero and behaves as 
Gi (A) = A~^, while all other observables are equal to 



zero within error range. It should be noted that the el- 
ement {n, ±} with n > 1 are generated by our random 
process, but their contributions to Gn at n > 1 cancel on 
average, and as a result G„ (A) = for n > 1. 

A similar stochastic strong-coupling solution can be 
easily constructed for other models with U (N) field vari- 
ables, for example, for the large- A'' U (N) non-Abelian 
lattice gauge theory. In this case the desired random 
process will be a minor modification of the stochastic 
solution of the Weingarten model However, such a 
strong-coupling solution is not very interesting, since the 
continuum physics should emerge in the weak-coupling 
limit. 



Weak-coupling regime 

At sufficiently small A the equations ^ cannot be 
transformed to the stochastic form by a simple rescal- 
ing of observables. In order to obtain a perturba- 
tive solution in this case, let us introduce the new 
Hermitian matrix variables X via the relation g — 
{1 + iaX) {1 ~ iaX) ^, where the value of a will be 
given later. First we have to express the Haar mea- 
sure on U (N) in terms oi X. It can be found from the 
invariant metric form Tr {dgdg~^^ = gAsdX^dX^ — 

4a2Tr ( (l + a^X^V' dX (l + a^X^V^ dx): 
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Now one can include the Haar measure (O into the action 
for the variables X and express integrals over g as inte- 
grals over X. The new action can be expanded in powers 
of X, thus leading to the Hermitian matrix model with 
an infinite number of terms in the potential: 
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Now we set a = A/ (2A -I- 8) so that the quadratic term in 
the action ([6]) becomes l/2TrAr^. It is interesting that 
another parameterizations of U (N) group manifold by 
Hermitian matrices (e.g. the exponential g — exp (iaX)) 
can in general lead to multi-trace terms in the action. 

The SD equations for the new action ^ involve the 
observables ^„ = (-^TrAT^"). Again, we rescale the 
variables as ^„ — J\fc"Wn to cast the SD equation in the 
stochastic form ([1]): 
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FIG. 2: Observables Gn (A) in the weak-coupling phase. 

Following [HI, one can devise a random process which 
samples the unknowns w„. The configuration space of 
the process is again the stack containing positive integers 
which label the observables ^„ and signs. At each time 
step one performs at random one of the following actions: 

Create : With probability (AAc) ^ push a new element 
{1, +} to the stack. 

Join : With probability TV pop the two elements 
{m, si} and {n, S2} and push a new element 
{m + n + 1, S1S2}. 

Increase : With probability 2/c increase the topmost 
element in the stack by one. 

Decrease : With probability 2c^a^+^ (4fc + 4 + A) de- 
crease the topmost element {n, ±} in the stack by 
n > /c > and change its sign if k is even. 

Restart : Otherwise restart with a stack containing a 
single element {1,-1-}. 

The total probability of all actions can be made less than 
one ii \ <\w ~ 0.2. 

In order to calculate the observables G„, one should ex- 
pand them in powers of X: Tr [g (A)" -I- g^" [X)) ~ 
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1 — Gn can be then found from 



k=l 



the averages {{Afc'^gn^ sign (fc) ) ) w.r.t. the above de- 
scribed random process. Here k is the topmost element 



in the stack and sign (k) is the corresponding sign. The 
observables Gn (A) are compared with the exact solution 
on Fig. [2] for A = 0.2 w A^, and for A = 0.1. In order 
to illustrate the efficiency of the algorithm, on Fig. [3] we 
plot the contributions of different k (positive and nega- 
tive separately) to 1 — Gn (A) at A = 0.2 « X^. Solid 
lines are the exact values G„ (A). It should be noted 
that since a ~ A, each such contribution corresponds 
to the term of order A*"' in the weak-coupling expansion 
of Gn (A). Gi (A) is saturated almost completely by the 
lowest-order term, while for Gio (A) large contributions 
of different k cancel to yield a numerically small value. 
Contributions of negative k are numerically small in both 
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FIG. 3: Contributions of different terms of the weak-couphng 
expansion to Gn [Xw), for lO" iterations of the random pro- 
cess. 



cases. Thus the strength of the "sign problem" for our 
random process actually depends on the choice of ob- 
servables. We also note that no critical slowing down 
near A^.^^ was observed neither in the weak- nor in the 
strong-coupling phases. 
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